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< : ABSTRACT 

CN . With the approaching start-up of the experiments at LHC, the urgency to quantify systematic 

uncertainties of the generators, used in the interpretation of the data, is becoming pressing. The 

>• i PHOTOS Monte Carlo program is often used for the simulation of experimental, selection- 

C^ I sensitive, QED radiative corrections in decays of Z bosons and other heavy resonances and 

CNl ' particles. Thanks to its complete phase-space coverage it is possible, with no approximations 

O 



for any decay channel, to implement the matrix-element. The present paper will be devoted to 
^ , those parts of the next-to-leading order corrections for Z decays which are normally missing in 

O I PHOTOS. The analytical form of the exact and truncated, standard, kernel used in PHOTOS 

^ I will be explicitly given. The correction, being the ratio of the exact to the approximate kernel, 

can be activated as an optional contribution to the internal weight of PHOTOS. 

To calculate the weight, the information on the effective Bom-level Z/y* couplings and even 
^ . directions of the incoming beams, is needed. A universal implementation would have made the 

PHOTOS solution less modular and less convenient for the users. That is why, for the time 

being, we will keep the correcting weight as an extra option, available for special tests only. 
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;h ' We will quantify the numerical effect of the approximation with the help of a multitude 

of distributions. The numerical size of the effect is in general below 0.1%; however, in some 
corners of the phase-space (well defined and contributing less than 0.5% to the total rate), it 
may reach up to about 20% of their relative size. 
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1. Introduction 



Analysing the data from high-energy physics experiments, we try to solve the "experiment 
= theory" equation. This non-trivial task requires many different effects to be considered simul- 
taneously. From the experimental side, these are mainly detector acceptance and cuts, which 
are dictated by the construction and physical properties of the detector: the shapes of distribu- 
tions may be distorted by, say, misidentification and residual background contamination; these 
effects need to be discriminated in an appropriate and well-controlled way. From the theoretical 
side, all effects of known physics have to be included in predictions as well. Only then can 
experimental data and theoretical predictions be confronted to determine numerical values of 
some coupling constants or effects of new physics (to be discovered). 

A well-defined class of theoretical effects contains the QED radiative corrections. PHOTOS 
is a universal Monte Carlo algorithm that simulates the effects of these corrections in decays of 
particles and resonances. It is a project with a rather long history: the first version was released 
in 1991 [1], followed by version 2.0 [2] (double emission, threshold terms for fermions). The 
package is in wide use [3] : it was applied as a precision simulation tool for W mass measurement 
at the Tevatron [4] and LEP [5,6], and for CKM matrix measurements in decays of K and B 
resonances (NA48 [7], KTeV [8], Belle [9], BaBar [10] and at Fermilab [11]). 

Throughout the years the core algorithm for the generation of 0{a) corrections did not 
change much; however, its precision, applicability to various processes, and numerical stability 
improved significantly. New features, such as multiple photon radiation or interference effects 
for all possible decays, were also introduced. 

Growing interest in the algorithm expressed by the experimental collaborations (including 
the future LHC experiments) was a motivation to perform a more detailed study of the potential 
and precision of the PHOTOS algorithm. The present paper is the third in the series [12, 13]. It 
is devoted to the Z boson decay and to simplifications in the matrix element used in PHOTOS 
for that channel. We also explore the limitations originating from compromises introduced into 
PHOTOS bremsstrahlung kernels, which assured convenience of use; no process-dependent 
weight need be involved. 

In that respect, the study of the PHOTOS matrix element can be understood as a part of 
the on-going effort to find the practical solutions of the improved expansions. Some aspects 
of our solution resemble those of classical exclusive exponentiation as described in [14, 15]; in 
an other, the parton shower may be identified. The solution may be understood as a rearrange- 
ment of the QED perturbation expansion, yet this point will not be discussed here. Instead, 
let us point to some similarities of the PHOTOS solution to the methods discussed elsewhere: 
interaction picture of Quantum Mechanics, expansion of special functions around asymptotic 
solutions [16] or in field theory; see eg. [17]. In PHOTOS the expansion is performed in terms 
of multidimensional operators. 

The paper is organized as follows. In Section 2 the main properties used in the PHOTOS 
design and, in particular, the analytical form of the (NLO) weight, necessary to introduce the 
complete first-order matrix-element, are presented. It is also explained there how the com- 
plete matrix elements break the requirement of separation of the calculation of the final-state 



bremsstrahlung from the properties of the Bom-level matrix elements and the phenomena af- 
fecting the Z production. To support the discussion and visualize the results, a multitude of 
numerical comparisons and tests will be presented. Section 3 provides the definition of the 
method used in those comparisons. The method is particularly suitable to visualize the results 
in the non-coUinear regions of the phase space. Section 4 presents numerical tests performed 
at fixed first order of the QED expansion. Since PHOTOS uses the same building block for a 
part of the single-photon generation algorithm and for the multiple bremsstrahlung, the results 
presented in this section have implications for the multiple-photon option of PHOTOS. Section 
5 addresses these aspects of the program construction which are relevant to the use of the NLO 
weight in the multiple-photon option. Section 6 collects the results of the tests performed for 
the programs run with multiple-photon emission. Summary, Section 7, closes the paper. 



2. Phase space and matrix element 



To discuss the implementation of the complete first-order QED radiative corrections in Z 
decay, we must start with the complete parametrization of the phase space. 

Let us start with the explicit expression for the parametrisation of the {n+ l)-body final 
state in the decay of an object of four-momentum P. To define iterative relations, let us denote 
the four-momenta of the first n decay products as kj, and the last {n + l)-th as q. For the case 
discussed here, the (n + l)-th particle will always be the real and massless photon. However, 
the parametrization does not rely on this assumption and, in principle, can be applied to define 
other formulas for the phase space, such as the emission of a (massive) pion, and could even be 
extended to the case of emission of pairs of heavy particles. In later steps of our construction 
the fact that photons are massless and the related properties of QED matrix elements will of 
course be used. 

In the following, the notations from refs. [18, 19] will be used. We will, however, not rely on 
any particular results of those papers and only point to the more detailed presentations of other, 
nonetheless quite similar, options for the exact n-body phase-space parametrization to the one 
presented here. 

Let us define the element of Lorentz-invariant phasespace (Lips) as follows: 
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Extra integration variables, the four- vector p, compensated with 6^ {p — £" k^ , are first intro- 
duced. In the next step, another integration variable Mi and 5(p^ ~^?) ^^e introduced. The 
element of the phase-space integration may thus be transformed into: 
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The part of the phase-space Jacobian corresponding to the integration over the direction and 
the energy of the last particle (or invariant mass of the remaining system) is explicitly given; 
X(a, b, c) = a^ + b^ + c^ — lab — lac — Ibc. The integration over the angles is defined in the 
rest frame of n + 1 particles; the integration over the invariant mass M\ is limited by the phase- 
space boundaries. There is no need to choose the axes with respect to which the angles are 
oriented; we will not elaborate on that point here, as details can be found in Ref. [1]. Formula 
(EJ) may be iterated to provide parametrization of the phase space with an arbitrary number of 
final-state particles. The question of the orientation of the reference frames used to define the 
angles and the order of the choice for limits in M, integrations, becomes particularly complex 
then; our choice is described in ref. [2] . Since nothing new was introduced for the purpose of the 
present study we will not discuss this interesting point further. Except for the details mentioned 
above, the choice we made for the phase-space organization is the same as in FOWL [20], 
TAUOLA [19], and probably many other generators. 

To simplify the formula for the phase space, let us finally take advantage of the zero mass 
of the photon. The invariant mass of the system of all particles but the first one may be replaced 
by the energy of the first one (defined in the P rest frame). The phase- space formula can then 
be written as: 
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If we had / photons accompanying n other particles, the factor in square brackets would be it- 
erated. A statistical factor jj would complete the formula for the phase-space parametrization, 
which is quite similar to the formal expansion of the exponent^ . The last formula, supplemented 

'The exact form of the functional exponent is achieved if the four- vector p is replaced by P in formula (|3}- In 
this way the tangent space for the (« + l)-body phase space can be constructed. We use that space, together with 
an eikonal-like form of the matrix element (emissions from individual final-state charged products are treated as 



with the definition of the orientation of the angles, is used to define the full kinematic configura- 
tion of the event. The four-momenta of all final-state particles may now be constructed from the 
angles and energies (fcy,) of the photons, and the angles and masses of the other decay products. 

Similarly, an inverse operation may be performed; the energies and angles for the parametriza- 
tion could be reconstructed from the four-vectors (even though the parametrization was not 
necessarily used in the previous generation steps). The phase-space Jacobians may be easily 
calculated as well. By replacing dLipSn{p —^ ki...k„) in formula ^ by dLips„{P — > k\...kn) 
we obtain a parametrization where the photons do not affect the construction of other particles' 
momenta. This operation could be considered as treating the photon in an approximation valid 
only in the soft photon limit. This, however, does not need to be the case. In the first step, the 
photon may be constructed with an arbitrarily large momentum, as nothing else depends on it. 
The kinematical variables of the photon are generated with the help of the distribution defined 
by the factor \kydkyd co?,Qd^j7Y^]^ ; which provides the photon variables of the tangent space. 
Fully constructed with four-momenta, an event of the n-body decay can be turned back into a 
representation of angles and invariant masses. In the final step these angular variables, together 
with those of the photon, can be used to define a new event in the (n + l)-body phase space. 
In the case when the new kinematical variables do not fit the limits of available {n + l)-body 
phase space the new event should be rejected and the original configuration (in the n-hoAy phase 
space) kept. An important property of the algorithm presented here is the full coverage of the 
{n + l)-body phase-space being assured. In this procedure, the difference between n-body and 
n-\-\ body phase-space Jacobians can be calculated in an unambiguous way and introduced in 
the same rejection step as for the phase-space limits^. 

The features and transformations of the phase- space parametrization presented here are at 
the heart of the construction of the PHOTOS kinematics and have been used since its begin- 
ning. To complete the generation of photons, the exact phase-space parametrization must be 
completed with a matrix element, with both virtual and real QED corrections included. Careful 
regularization of soft singularities must be performed^. 

independent), for the construction of the crude distribution of photon emission probabiHty. Note that, in this space, 
the photons' four-momenta are unconstrained by energy-momentum conservation. The limits on the energies of 
the photons are arbitrary. We checked that, at the lO^"* precision level, the results obtained from our simulations 
do not depend on the particular choice. We leave the underlying formal aspect of the algorithm to future papers. 

^The effects of matrix elements, including those of virtual corrections, have to be introduced at this stage as 
well. They are indispensable, for example, to calculate relative probabilities of configurations with distinct number 
of final-state particles. 

^Volumes of the partial width attributed to the configurations with «, « + 1 particles, etc., have to be normalized 
to the total width, both at the level of the tangent and the correct (final) phase space. 
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In the standard version of PHOTOS, as published in [1,2], the expression 
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where: 0+ = Z(p+,^+), 0_ = Z(p_,^_), 

0Y= Z(y,;U^) is defined in (/j+,/j^)-pair rest frame. 
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is used for the real-photon matrix element. The virtual corrections are requested to be such 
that the total decay rate remains unchanged after complete QED corrections are included. The 
expression, without approximation, reads: 
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The combined effect of the virtual and real corrections on the total rate is its increase by a factor 
The notation from ref. [21] are used: 



s = 2p+-p- 



2q+ ■ q^ 



t = 2p+-q+, t' = 2p+-q- 
u = 2/?+ ■ q , u' = 2-- q^ 

k±=q±- k, Xk = 2E^ 



(6) 

This paper collects complete first-order radiative corrections for the process e'^e~ -^ iJ~^iJ~{y). 
Final-state bremsstrahlung constitutes part of these results, where nonetheless reference is made 
to the incoming electron beam momenta. 

The A term encapsulates final-state mass -dependent terms, /?+, p-, q+, q^, k denote four- 
momenta of: incoming e+, e^ , outcoming /j+, /j and the bremsstrahlung photon respectively. 
Expression Q is explicitly taken from ref. [21], on Monte Carlo MUSTRAAL, which is where 
the interested reader will find the details of the definitions of variables and expressions such as 
A, used in the formulae @ and ©. 

The ratio of Q to @ constitutes the basic element of upgrading PHOTOS functionality to 
the complete first order^. Nothing needs to be changed in the phase-space parametrization. The 

'*This is only true for PHOTOS being run at first order. For the multiple-photon radiation option, the iteration 
of the single-photon emission kernel (and thus also its weight) is performed; see section 5. 



effects of the virtual corrections have to be included as well and have to be properly introduced 
in the normalization. The expression for the correcting weight could be chosen^ simply as 
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For the purpose of constructing a Monte Carlo algorithm, however, it is more convenient to sep- 
arate it into a sum of two generation branches (with slightly different angular variable mapping). 
Then, the expression for the distribution and those for the weight take the form 
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At this point, let us make the following remark. Event though the introduction of the NLO 
weight into PHOTOS is trivial, the developed approximation [1] at the heart of PHOTOS design 
is not. It enabled universality of the program^. Simplification was not necessary to attribute the 
generation of bremsstrahlung photons to individual charged particles^. The separation holds for 
the complete NLO as well. The simplified emission kernel, which we used for other decays as 
well, reads: 
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^Alternatively, a factor (1 + | ^) can be included in the definition of the crude distribution. 

^Indeed, after inspection, the differences between formulae and (|5jl are quite significant. The exact ex- 
pression does not allow a transfer of the complete Born-level angular dependence to the host generator In the 
correction weight, the two contributions, one depending on the angle 0+ and another on 0_, have to be simultane- 
ously included. The dependence on the Born-level (effective) couplings thus need to be known at the level of the 
calculation of the final-state bremsstrahlung weight. This would make the modular structure of PHOTOS design 
more difficult to keep. Also, the direction of the (effective) beam need to be provided for the calculation of 0+ and 
0- angles. This exhibits another difficulty in the separation of the final-state bremsstrahlung and the dynamics of 
the initial state for Z/y* production. 

^For other decays, it will probably not be necessary to find an explicit form of such NLO separation. Starting 
from the NNLO, such separation was shown to be impossible [22] anyway. 



It depends on the spin and charge of the "emitting particle" only^. It does not depend on the 
properties of the other decay products, which only define the phase- space limits. To obtain the 
universal form of the photon emission kernel, the interference was eliminated with the help of 
the factor, either n , a r^ \ or /, ^ ^ ^ . The interference is recovered later, using the weight 

' (l + pCOSfc)y) (1 — pCOSBy) 

given in formula (17) of ref. [2], that is with approximation. On the other hand, having paid the 
price of the approximated solution, both the kernel and the interference weight can then be used 
for decay of any particle or resonance. 

In our present study the analytical expression for the matrix element for the e+e^ -^ Z'^/'f — > 
fj^fj~ is used (and compared to the approximated, yet process-independent, solution of the stan- 
dard PHOTOS). If necessary, a matrix element for other decay processes (if available) could be 
used. 



3. Method used in numerical tests 



In the comparison of the multitude of final states generated, at different levels of physics 
sophistication by two distinct Monte Carlo programs, the choice of a method is of great impor- 
tance. To compare the Monte Carlo programs it is quite common to present the distributions 
generated by the programs superimposed on a single plot, often in logarithmic scale. Such a 
method was used, for example, in ref. [23]. The method is unquestionably sufficient, if one's 
interest is limited to, say, the coUinear content of the results or other distributions of the intrin- 
sically logarithmic type. 

For instance, applying this method for the comparison of the total energy carried out by 
all bremsstrahlung photons, we would obtain a distribution such as those presented in Fig. [l] 
We could conclude that there is excellent agreement, and the non-leading effects, which are 
essential for estimating systematic errors for generators like PHOTOS, would be marginalized 
in the presentation. This is also the case if for the same distribution (see right-hand side of Fig. 
[B the method [13] based on MC-TESTER is used. The distributions are indeed dominated by 
the coUinear content of the programs! For other distributions, sensitive to second-order matrix- 
element parts, missing in PHOTOS, the differences would become visible on the plots obtained 
normally from MC-TESTER. That is why, in the present paper, we will keep to that class of 
comparison plots. The comparisons are automated and standardized. This not only reduces the 
time needed for debugging the tests, but also allows for easy cross-comparisons of the results 
presented in our consecutive papers. 

For a selected decay process, such as the Z/y* decay, the four-momenta of the decay prod- 
ucts and their flavours are extracted from the event record in an automated way (thereby limiting 
the effort of setting up the appropriate analysis code and also the risk of accidental errors). The 

**In the original PHOTOS documentation we called this well-controlled truncation of the kernel a "property 
(such) that leading-log (coUinear) and infrared limits are properly reproduced". This explanation turned out to be 
misleading for many readers. One can get biased, and expect the collinear approximation not only for the kernel, 
but for the whole design of the algorithm. This would be a serious limitation of our program design if indeed, 
as suggested in ref. [23], "PHOTOS was based on collinear approximation". Fortunately it is not the case. Such 
confusion was not a concern for the users, until now: precision requirements were not as high. 



Figure 1: A typical plot for comparisons, as described in ref. [23]. We use it to illustrate our method. 
The histograms presented in the left and right plots have identical content, and show a comparison of 
KKMC [14] and PHOTOS used in Z decay. The total energy carried out by all final-state photons is 
presented. Red (darker grey) colour represents the results of KKMC, green (lighter grey) of PHOTOS 
with the NLO weight activated. Samples of 10' events were used in this comparison. The results are 
overwhelmed by the collinear/soft content of the predictions. 

If the W instead ofZ decay was chosen, it would not be the case. There, NLO effects would be dominant 
for the part of the spectrum above My^ jl. However, in that case, we would not profit from the second- 
order matrix element Monte Carlo, available for tests. 





(a) A logarithmic scale is used. Excellent agree- 
ment between the two programs is visible all over 
the energy range from to Mz- The presence of 
the two lines can be spotted at the high end of the 
spectrum, mainly thanks to the statistical errors. 
The kink at the limit of the phase-space, where sin- 
gle hard-photon configuration ceases to contribute, 
dominates the content of the whole picture. 



(b) The method of MC-TESTER is applied. Both in- 
dividual distributions, from PHOTOS and KKMC, 
are presented, but overlap. The samples populate 
the first few bins of the histograms. The differences 
would normally be visible on the black histogram, 
which presents the ratio of the results from PHO- 
TOS and KKMC. The agreement is perfect all over 
the spectrum. No structure can be spotted in the 
vicinity of the kink (total photon energy, close to 
half of the Z mass). If present, the structure of pos- 
sible differences, would be well separated from this 
of the shapes themselves. 



decay events obtained that way are classified in distinct decay channels, according to the par- 
ticles present in the final state. The histograms of all possible invariant masses, which can be 
formed from the decay products, are defined and filled for each identified decay channel. At 
the end of the run they are stored in output files. Two output files (from distinct runs of event 
generators instrumented with MC-TESTER) are then analysed, and the results are presented in 
a form visualized as a "booklet" made of plots and summary tables. The user is given some 
general information concerning the comparison of the two runs with different Monte Carlo gen- 
erators, a list of the decay channels with their branching fractions, and the maximum values (for 
each decay channel) of the shape difference parameter (SDP)^. 

For each decay channel the plots of histogrammed values are then included; each plot 
presents two distributions from the two distinctive runs and a curve, which is the ratio of the 
two normalized distributions. The value of the SDP is also printed for each plot. In practice, 
as in paper [13], the histograms obtained from the compared programs will often overlap. The 
differences will then be visible only in the plot of the ratio of histograms. 

The testing approach implemented in MC-TESTER could be used directly in the case of 
validation of the TAUOLA package. Nevertheless, for the purpose of studies presented here, it 
needed an extension. It is necessary, according to the particular method of handling soft photon 
cancellations, to consistently treat the soft final-state QED bremsstrahlung photons, which may 
or may not be present in the event. If results of different programs were compared blindly, 
ambiguities due to differences in the treatment of the soft emission region and of the different 
boundaries for the photon phase-space (integrated analytically) would arise. To prevent these 
ambiguities, the most convenient solution was to introduce a technical regulator in the test itself. 

For our comparisons to make physical sense and remain automatic, we had to remove the 
softest photons from the final states. We defined zero-, one-, and two-photon topologies in 
the following way: we called the event "zero photon" if there was no photon of energy (in a 
decaying particle's rest frame) larger than ^test- The "one-photon" event had to have one (and 
only one) photon of energy larger than ^test- If there were more than one such photons, we 
called it a "two-photon" event. In the case where there were more than two photons of energy 
larger than ^test^ we considered only the two most energetic ones, and treated the remaining 
(softer) ones as if they had not passed the fitest threshold. For all the photons that did not pass 
the ^test threshold we summed their four- momenta with the momentum of the outgoing fermion 
of smaller angular separation. With the help of our test we divide the phase space for two 
fermions and an arbitrary number of photons into slots of 0, 1 and 2 distinguished final-state 
photons. 

In the paper we will use two variants of this test definition: testl and test!. The test! is 
exactly as explained above. In testl, only one photon (the most energetic one) will be accepted. 
The free parameter, Etest is chosen to be 1 GeV for all results presented in this paper. 

Systematic histogramming of all possible invariant masses that can be constructed from a 



'^The shape difference parameter, defined in [24], quantifies the difference in shape of the histograms coming 
from the two runs being compared. The SDP value is calculated separately for each histogrammed mass: it 
quantifies the exclusive surface between the (normalized to unity) corresponding histograms obtained from the two 
runs. The effects of statistical fluctuations are appropriately subtracted. The maximum SDP over all distributions 
for a given decay channel is taken and printed in the table. 



combination of final-state four- vectors, and storing them as one-dimensional histograms, does 
not define a test of ultimate sensitivity. The method is blind to the P-parity- sensitive effects, 
important for x lepton physics, for instance. Also important effects, such as coherence between 
the photons are to a large degree washed out. Nonetheless we believe that the advantages of the 
method are prevailing, and we decided to use it in this study. 



4. Results of the tests performed at first order 



Let us start with a comparison of PHOTOS and KORALZ [25], both run at first order and 
without exponentiation. For KORALZ, the complete first-order matrix element, as in Section 
2, is used. The results from KORALZ are given by a red (darker grey) line and from PHOTOS 
by a green (lighter grey) line. (In the presentation of the results we use the colour coding 
consistently in the plots and the summary tables, following the methodology of MC-TESTER). 
The lines overlap almost completely on all plots and only the ratios of the distributions shown 
as black histograms indicate that there is some difference. The actual plots for MC-TESTER 
comparison, are prepended with a summary table giving the fractions of event with and without 
photons (of energy above 1 GeV). In all comparisons samples of 10^ events were used. 



Decay channel 


Branching ratio ± rough errors 
KORALZ PHOTOS 


Max. SDP 


Z^-AZ-A/^ 


82.5137 ± 0.0091% 82.3622 ± 0.0091% 


0.00000 


Z^^^-^+y 


17.4863 ± 0.0042% 17.6378 ± 0.0042% 


0.00534 



As can be seen, the difference in the fraction of events with photon (of energy above 1 GeV) 
is about 0. 15%. Although noticeable (thanks to our method), this is not a large discrepancy. Let 
us now turn to the distributions. 



Decay Channel: Z^ —^ ju jj^y 
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Comparison of Mass(1 ) of mu- mu-f gamma in channel ZO => mu- mu+ gamma 
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Analysing the values of the SDP, printed in the upper right comers of the plots, we conclude 
that the surfaces between green and red histograms (both normalized to unity before calculation 
of the ratio) differ by at most 0.005; this low value quantifies the fact that the histograms for the 
two programs overlap almost completely. The ratio of the two lines (black histogram) nonethe- 
less reveals the difference, which is located at the far end of the spectra, sparsely populated 
by configurations with photons of extremely large energies and away from the direction of the 
muons. 

Even though the agreement is amazing, it is yet improved once the NLO term is included 
into the correction weight of PHOTOS: the differences disappear below the statistical error of 
10^ event samples and they are not noticeable in the histograms' ratio curve either! In the 
following table and figures, we collect the results as previously discussed, but for the runs with 
the NLO correcting weight activated in PHOTOS. 



Decay channel 


Branching ratio ± rough errors 
KORALZ PHOTOS 


Max. SDP 


Z^^A/V 


82.5110 ± 0.0091% 82.5074 ± 0.0091% 


0.00000 


Z" -^ A'+A'^Y 


17.4890 ± 0.0042% 17.4926 ± 0.0042% 


0.00000 
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Decay Channel: iP -^ ju^ju y 



Comparison of Mas5(1) of mu-f mu- in channel ZO => gamma mu-f mu- I ^Qp 
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The agreement for the branching fractions (of events with and without photons of energy 
larger than 1 GeV) is better than 0.01% now! This excellent agreement indeed confirms that the 
theoretical effects missing in the standard version of PHOTOS are negligibly small. It is equally 
important that it provides, a powerful technical test of the generator. The kinematical variables 
used in PHOTOS differ from those of KORALZ; four-vectors are used instead of angles to 
parametrize the intermediate steps of the generation. The differences could have indicated, 
say, consequences of aggregation of rounding errors. Keeping in mind that similar levels of 
agreement for muons was achieved for the multiphoton version of PHOTOS and KKMC in the 
case of Z ^ e+e^ decay, we can confidently claim that PHOTOS has numerical stability under 
control. This was not the case for the early versions of the program, and reaching that level of 
technical reliability required a major effort. 

5. Algorithm for multiple-photon generation 

Before presentation of the results for multiple-photon generation from PHOTOS, let us com- 
ment on those technical details of the PHOTOS algorithm, that are important in the implemen- 
tation of the NLO contribution to the correcting weight. The iteration algorithm, as explained 
first in ref. [2], and recently also in refs. [12, 13], did not require changes for the case of mul- 
tiple photon generation. Nevertheless, the following details have to be clarified for the proper 
implementation of the NLO weight, given by formula (El). All identical terms present in the 
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numerator and denominator, expressions ^ and (El) respectively, were cancelled out at the ana- 
lytical level^*^. The weight © is always calculated for the single-photon configuration. If there 
are other photons generated in the previous steps of the iteration, their momenta are absorbed 
into the momenta of the final-state fermions. The constraint on the direction and the opening 
angle between the photon under consideration and the direction of the charged emitter [1], is 
assured. 



6. Numerical results of the tests performed with multiple-photon radiation 

Let us now turn to the tests of PHOTOS running in multiple-photon option. For that purpose 
we will use testl as defined in Section 3, and samples of 10^ events generated from KKMC with 
exponentiation and the second-order matrix element and the multiple-photon radiation version 
of PHOTOS, without NLO terms. The results are included in the table and plots below. 



Decay channel 


Branching ratio ± rough errors 
KKMC PHOTOS 


Max. SDP 


Z'^iu-iu^ 


83.9176 ± 0.0092% 83.8372 ± 0.0092% 


0.00000 


Z^^^-^+y 


16.0824 ± 0.0040% 16.1628 ± 0.0040% 


0.00409 



Decay Channel: Z^ — > ^ jn^y 



Comparison of Ma5s(1) of mu- mu-t- in channel ZO => mu- mu+ gamma I (>l])p 

0.00409 



Comparison of Mass{1) of oiu- gamma in channel ZO = 




gamma | gpp 

0.0025 




20 40 



'"The cancelled-out terms could have been calculated using slightly different kinematical variables; the differ- 
ences would have appeared only in the case of more than one hard photon present in the final state. In such a case 
the ratio of the terms would not be equal to 1 . These effects generally go beyond the NLO, and in fact our choice 
was motivated by the comparisons with the second-order matrix-element calculation, but without necessary details. 
That is why an appropriate discussion of this choice would require detailed presentation of the second-order matrix 
element. It would have to be similar, for example, to the discussion of the extrapolation procedure as that described 
inref. [26]. 
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Comparison of Mass(1) ot mLj+ gamma in channel ZG = 




Comparison of Mass(1 ) of mu- mu-f gamma in cfiannel ZO =; 
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The pattern of differences between the results of PHOTOS and KKMC runs resembles the 
one present in the plots for the comparisons performed at the first order. Again, the black 
curves of the histogram ratios for the KKMC and PHOTOS results are not consistent with 1 
for configurations with hard photons and in the regions where histograms are nearly at 0. The 
differences reach few per cent in the comers of the phase space contributing few per mille to 
the total rate. The discrepancies are again smaller than 0.1% with respect to the total rate. 

Once the NLO weight in PHOTOS is activated, the already small differences become even 
smaller, by a factor of about 50, measured with the SDP The differences are practically with 
10^ samples, which can be seen in the table and the plots below. 



Decay channel 


Branching ratio ± rough errors 
KKMC PHOTOS 


Max. SDP 


Z"^A'"A'+ 


83.9176 ± 0.0092% 83.9312 ± 0.0092% 


0.00000 


Z^^lj-lj+y 


16.0824 ± 0.0040% 16.0688 ± 0.0040% 


0.00003 



Decay Channel: Z' 







jj jj+y 



Comparison of Ma5s(1) of mu- mu-t- in cfiannel ZO => mu- mu-t- gamma I ^J)p 
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Comparison of Mass(1 ) of mu- mu-f gamma in channel ZO =; 



JsDP 



ixlO 

-10 

^8 

^6 

^4 

-2 



20 40 60 



80 100 



This confirms tliat tlie main source of residual discrepancy between KKMC and PHOTOS, 
both running in exponentiated versions, was due to the NLO term missing in the previous group 
of results and activated now. 

Let us now turn to test!, where configurations of up to 2 hard photons are analysed. To 
this end, we will present at first the results of a comparison of the standard multiple-radiation 
version of PHOTOS with that of KKMC, followed by the comparison of multiple-radiation 
version PHOTOS with NLO weight and KKMC (again with second-order matrix element and 
exponentiation). 



Decay channel 


Branching ratio ± rough errors 
KKMC PHOTOS 


Max. SDP 


Z"-A'-A'+ 


83.9177 ± 0.0092% 83.8372 ± 0.0092% 


0.00000 


IP -^ jj-jj+y 


14.8164 ± 0.0038% 14.8676 ± 0.0039% 


0.00232 


Z^^/j-jj+YY 


1.2659 ±0.0011% 1.2952 ±0.0011% 


0.00918 



Decay Channel: Z^ ^ ju ju^y 



Comparison of Ma5s(1) of mu- mu-t- in channel ZO => mu- mu-t- gamma I (>1J)P 

0.00232 
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Comparison of Mass(1) ot mLj+ gamma in channel ZG = 
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Decay Channel: Z^ ^ fi jn^yy 




Comparison of Mass(1 ) of mu- gamma in channel ZO =: 
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One can see that already for the standard PHOTOS the agreement is good. Residual defi- 
ciencies are small for both slots of the phase space: single hard photon, and two hard-photon. 
Let us present now what kind of changes the inclusion of NLO terms in PHOTOS brings to the 
results of test2. 



Decay channel 


Branching eatio ± rough errors 
KKMC PHOTOS 


Max. SDP 


Z"-A'-A'+ 


83.9177 ± 0.0092% 83.9303 ± 0.0092% 


0.00000 


Z" -^ p-p+y 


14.8164 ± 0.0038% 14.7829 ± 0.0038% 


0.00005 


Z^^/j-lj+YY 


1.2659 ±0.0011% 1.2868 ±0.0011% 


0.00293 
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Decay Channel: 7? ^ \i ^+y 




Comparison of Mass(1) of mu' gamma in channel ZO = 
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Decay Channel: Z^ 
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Comparison of Mass(1 ) ot mu- gar 
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For the single hard-photon distributions, the differences diminished significantly, again a 
factor of about 50! This was to be expected. Even for the distributions of the phase-space 
slot with two hard photons, the differences diminished. The SDP decreased by a factor of 
about 3. This is not as striking as for the single hard-photon configuration, but it is of no 
surprise: the complete second-order matrix element is missing. The improvement by a factor of 
3 provides, however, a strong indication that the algorithm of iteration used in the generation of 
consecutive photons work well from the point of view of NNLL level as well. The acoplanarity 
plots presented in the previous paper [13] also demonstrated some of the NNLL aspect of the 
algorithm. That is why we are not going to discuss this point here, but we would rather leave it 
to future discussion of the NNLL content of our algorithm: this aspect goes beyond the purpose 
of the present paper and the scope of interest of most of PHOTOS users. 

Much as is described in the present paper, a new contribution to the PHOTOS correcting 
weight would be needed for the NNLL case. However no changes in the phase-space algorithm 
would a priori be required. The techniques of gauge-invariant separation of the amplitudes 
into parts, as used for instance in ref. [26], will probably be necessary. They proved to be 
instrumental in the implementation of the second-order matrix elements for e+e -^ VgVeYY 
into KKMC. The exclusive exponentiation scheme of the KKMC Monte Carlo is prepared for 
j'-channel processes. 

We have to admit that once the NLO terms are switched on in PHOTOS, the difference 
between its results and those of the second-order matrix-element generator KKMC are at the 
limit of being recognized, even if samples of 10^ events are used. For the case of the two-photon 
test, differences due to the missing second-order matrix element in PHOTOS can be observed; 
yet they are too small, and don't have enough structure, to understand their possible origin. 
A part of the differences may even originate from the third-order LL terms (after integration), 
which are missing in KKMC but generated in PHOTOS in the process of iteration. 

In any case, even for PHOTOS running with standard options, the differences affect only 
a tiny fraction of the Z decay phase space. Thus, we do not consider it to be of much interest 
to continue the discussion of the missing terms. Nevertheless, from a more fundamental side, 
we are disappointed by the fact that the comparisons did not provide numerical insight into the 
structure of the differences. The particularly interesting aspects of the study in the context of 
the extension of the algorithm for QCD did not bring any constructive indications so far. 

7. Summary 



To quantify the size of the NLL effects, which are normally missing in PHOTOS, we rein- 
stalled them back into the program, using the original complete first-order expression for Z 
decay. After the NLO correcting weight was installed, the differences between PHOTOS and 
KORALZ were below the statistical error of 10^ events and for all the distributions used in 
the tests. Both PHOTOS and KORALZ were run at fixed first order without exponentiation. 
The agreement provided a technical cross-check test for the two simulations. For the case of 
multiple-photon radiation in PHOTOS, a comparison with the KKMC generator [14] (exponen- 
tiation and second-order matrix element used) was performed. The implementation of the NLO 
terms in PHOTOS indicated, in the results of our universal test, an improvement by a factor of 
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about 50 for the observables sensitive to a single hard photon in the final states and remained at 
the level of better than 0.1% on the total rate for all other cases we examined. Because of the 
smallness of the residual differences, it was difficult to understand their structure and origin in 
the final states with two hard photons. 

The improvement in the agreement due to the introduction of the NLO correcting weight 
came at a price. Even though the weight is analytically simple and generation of weight 1 events 
remained possible, the calculation of the weight required information on Born-level coupling 
constants of the intermediate Z/y*. Also, the direction of the beam was necessary in the calcu- 
lation of the weight. These requirements threatened the modular organization of the PHOTOS 
solution, as used in the large Monte Carlo generation chains of experimental collaborations. 
Numerically, the introduced improvements are rather small and the deficiencies of standard 
PHOTOS are localized in the comers of bremsstrahlung phase- space populated by photons of 
very high energies and angularily, well separated from the final-state muons. Those regions of 
the phase space weigh less than 0.005 to the total rate and the differences in that region approach 
20% of their size. The effects are thus less than 0.1% of the total rate of the Z decay to muons. 
That is why we do not think it justified to complicate the PHOTOS algorithm and to enable the 
use of the NLO correcting weight in the general case. 

The analysis presented here concentrates not only on the numerical results for the final- 
state bremsstrahlung in Z decay, but also on various aspects of a mathematical organization 
of the program for calculation of radiative corrections in Z production and decay. Separation 
of radiative corrections into parts: (i) embodied in effective couplings of the hard-scattering 
process, (ii) final-state QED bremsstrahlung, and (iii) initial-state bremsstrahlung, eventually 
with initial-state hadronic interactions, were mentioned as well. The effects of QED initial- 
final-state bremsstrahlung interference were to a large degree neglected. Such an approach is 
reasonable in the leading-pole approximation for the Z, but at a certain precision level the effects 
may need to be taken care of. For the time being the results of ref. [27] can be used instead. 

Thanks to the analytic form of the kernel used in PHOTOS for the single-photon emission, 
the analysis presented here may easily be extended to other decay channels, if high precision 
is required and a calculation of matrix element is available. The study for the case of B meson 
decay into a pair of 7r^(^^) is near completion [28]. In this case the questions of reliability of 
scalar QED for the calculation of photons of high pr, with respect to charged scalars, need to 
be addressed. A natural extension of the study of the systematic error in PHOTOS simulations, 
as presented here, would be the discussion of bremsstrahlung in W and Higgs boson decays. 

The decays of the W and Higgs bosons are probably the only ones where formal studies of 
the NL terms, similar to the ones presented in this paper, can be performed. As was the case 
with Z decays, those cases will also be limited to the leading-pole approximation. 

For other decay channels, the correction weight can be applied as well; however, in most 
cases the part of the weight going beyond soft and/or coUinear regimes may need to be con- 
structed with the help of the fits to the data. Let us stress that the unique design of PHOTOS, 
enabling the use of the same kernel for multiple-photon radiation (exponentiated) mode and at 
fixed first, second, third, and fourth orders, establishes a convenient environment for such fits of 
form factors to the data. At the same time the analytical form of these form factors can be taken 
from the first-order analytic calculations based on effective theory, or from any other model. 
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On the technical level it is worth mentioning that the NLO correcting weight of PHOTOS 
can be used as an internal correcting weight. 

Finally, let us stress that the approximations introduced in PHOTOS affect the matrix el- 
ements and not the phase space. The generation of the latter is based on the tangent space 
constructed from eikonal approximation but used also for hard photons, even of energies above 
the available maximum enforced by energy-momentum conservation. Only in the second step 
are phase-space constraints enforced. This is similar to the case of the classical exclusive expo- 
nentiation. The energy momentum constraints are introduced for each individual photon, step 
by step, and conformal symmetry is not used in that procedure. 

In principle, if necessary, complete higher-order matrix elements can be incorporated with 
the help of correcting weights. This point goes beyond the scope of the present paper. This is 
equally true for the possible extensions to simulations in QCD. 
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